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Quantum Monte Carlo (QMC) techniques are used to provide an approximation-free investigation 
of the phases of the one-dimensional attractive Hubbard Hamiltonian in the presence of population 
imbalance. The temperature at which the "Fulde-Ferrell-Larkin-Ovchinnikov" (FFLO) phase is 
destroyed by thermal fluctuations is determined as a function of the polarization. It is shown that 
the presence of a confining potential does not dramatically alter the FFLO regime, and that recent 
experiments on trapped atomic gases likely lie just within the stable temperature range. 
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I. INTRODUCTION 

Pair formation between fermions is one of the most rich 
and active fields of investigation in condensed matter sys- 
tems, occuring in contexts as varied as exciton formation 
in quantum well structures [l| to its most well-studied 
realization, the superconducting state 0. Indeed, the 
possibility of pair formation has many even wider im- 
plications, including the early speculation that it might 
dominate the nature of interior of neutron stars Q . By 
far the most commonly studied situation is one in which 
the populations of the two fermion species (normally up 
and down spins) are balanced, and Cooper pairs form 
with zero center of mass momentum Soon after the 
development of the BCS theory Q of superconductivity, 
the question of pair formation in polarized superconduct- 
ing systems, i.e. when the populations of the two spin 
states are imbalanced, was addressed independently b 
Fulde and Ferrel J5[ (FF) , Larkin and Ovchinnikov 
(LO) and Sarma |7J. The question was motivated by in- 
terest in the nature of superconductivity in the presence 
of a magnetic field, which induces spin polarization. 

FF and LO (henceforth FFLO) proposed similar, yet 
not identical, mechanisms whereby the Cooper pairs form 
with non-zero center of mass momentum equal to the dif- 
ference of the now-unequal Fermi momenta of the two 
species. As a consequence, in the FFLO scenario the or- 
der parameter is not homogeneous; the system has pair- 
rich regions separated by regions depleted in pairs but 
with an excess of the majority, unpaired fermion species. 
In the competing proposal by Sarma, the Fermi surfaces 
deform to allow pair formation with zero center of mass 
momentum. This means that the system remains uni- 
form; it is a homogeneous mixture of pairs and unpaired 
fermions from the majority population. In other words, 
in the FFLO mechanism, the momentum distribution of 
pairs has its peak at a momentum equal to the difference 
between the two Fermi momenta, fc pC ak = |&fi — kp2\- 
On the other hand, the Sarma mechanism would result 
in a pair momentum distribution with fc poa k = 0. Decid- 



ing between these two scenarios experimentally proved 
to be very difficult; the observation of the FFLO phase 
in solids was only achieved relatively recently in heavy 
fermion systems [8(. 

This question has recently taken on added importance. 
Interest in the astrophysical community has continued 
to grow and broaden. It is now believed that at ex- 
treme conditions of pressure, for example in the interior 
of supermassive stars, quark matter forms and pairing 
between the quarks could lead to color superconductivity 
More immediate experimental systems where these 
effects may be observed, and are indeed sought, are ultra- 
cold atoms confined in traps where two hyperfine states 
of fermionic atoms play the role of up and down spins. 
Such experiments have now reported the presence of pair- 
ing in the case of unequal populations [l(| EH in three- 
dimensional cigar shaped traps and in one dimensional 
traps [l2l] . However, the precise nature of the pairing has 
not yet been elucidated experimentally. 

On the theoretical side, much effort has been put into 
understanding the pairing mechanism. Calculations us- 
ing mean field theory [f3l423j . effective Lagrangian [24| 
and Bethe ansatz [25[ studies have been performed for 
the uniform system, with extensions to the trapped sys- 
tem using the local density approximation (LDA). Ex- 
tensive numerical work for the one dimensional system 
using Quantum Monte Carlo (QMC) HI HI, and the 
Density Matrix Renormalization Group (DMRG) [28l — 
Hf|) has demonstrated that, in the ground state, pop- 
ulation imbalance leads to a robust FFLO phase over a 
very wide range of polarization and interaction strengths. 

In addition to the above case of population imbalance, 
the case of mass imbalance has been addressed theoreti- 
cally [H, H(| and numerically [37j • This case is relevant to 
color superconductivity, where the quark masses are not 
equal, and also to ultra-cold 40 K- 6 Li atomic mixtures. 

The stability of the FFLO phase at finite temperatures 
remains an open question. Mean field calculations (32I — 
[H[ possibly shed some light but can be less reliable in low 
dimension where quantum fluctuations are large. This 
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question is of paramount importance for experiments, 
especially in trapped atomic systems, since difficulties 
in cooling fcrmionic atoms raise concerns about whether 
the currently attainable temperatures are low enough for 
a thorough investigation of FFLO physics. The present 
paper reports on QMC studies of this issue which provide 
an exact treatment of interaction effects on lattices of fi- 
nite size. The key result is that, when the polarization is 
sufficiently large, the paired phase can exist up to tem- 
peratures of order one tenth of the Fermi energy, even 
in the presence of a confining potential. Current exper- 
iments on trapped atoms are likely in this temperature 
range. 

The paper is organized as follows. In the next section 
we will discuss the model and the numerical methods 
used. In section III, we discuss our results for the phase 
diagram of the uniform one-dimensional system at finite 
temperature. In section IV we discuss the trapped one 
dimensional system followed by our conclusions in section 
V. 



II. MODEL AND METHODS 

In order to study the pairing mechanism of fermions 
in an optical lattice, we consider the one-dimensional 
fermionic Hubbard Hamiltonian, 
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where c\ a and c- are fermion creation and annihilation 
operators on lattice site i satisfying the usual anticom- 
mutation relation, {c ia ,Cj rT ,} — SijS^t. The fermionic 

species are labeled by a = 1,2 and n, CT = cl a c ig . is the 
corresponding number operator. The energy scale is set 
by taking the hopping parameter t = 1. The contact 
interaction strength is U which is negative since we are 
interested in pair formation in the attractive model. In 
the grand canonical ensemble, the populations are fixed 
by tuning the chemical potentials for the two species, \x\ 
and (i 2 , whereas in the canonical ensemble this term is 
absent and the populations are fixed simply by choos- 
ing the number of particles. As we discuss below, both 
ensembles will be used in this work. The last term de- 
scribes the confining harmonic trap which is centered at 
the midpoint, L/2, of the i-site lattice. We take periodic 
boundary conditions. 

Our Quantum Monte Carlo (QMC) results are ob- 
tained using two different methods: the Determinant 
QMC algorithm [H| (DQMC)and the Stochastic Green 
Function (SGF) technique jSl, HJ. In DQMC, the 



interaction term in the Hamiltonian, which is quar- 
tic in the fermionic operators, is decoupled via the 
Hubbard-Stratonovich (HS) transformation. This allows 
the fermion operators to be traced out, giving an ex- 
pression for the partition function as an integral over the 
configurations of the auxiliary HS field. We employ the 
grand canonical formulation of DQMC. The statistical 
weight of these configurations is given by the product of 
the spin up and down determinants, each of which can 
change sign depending on the HS configuration. When 
the populations are balanced, p\ = p 2 , it is well known 
that DQMC does not suffer from the sign problem for an 
attractive interaction (U < 0): The two determinants are 
identical and thus their product is always positive. This 
is no longer the case when the populations are imbal- 
anced, p\ ^ p 2 (or for U > 0), resulting in the appear- 
ance of the sign problem. The severity of this problem 
depends on several factors such as the size of the system, 
the coupling constant, the temperature and the imbal- 
ance in the populations. For the systems and parameters 
we consider in this work (see below) the average sign 
did not fall below 0.5, which allowed us to obtain good 
statistical precision in our simulations. We used this al- 
gorithm for the uniform system; it has the advantage of 
rather fast convergence rates, but to perform simulations 
at fixed particle numbers, the chemical potentials need 
to be tuned. A typical simulation for L = 50 sites and 
/3 = 30 with an imaginary time step At = 0.1 runs for 
about 4 to 5 days on a desktop computer and yields er- 
ror bars of the order of 0.5% for the peak of the pair 
momentum distribution, a central quantity of interest in 
this work. 

The SGF algorithm [3!| with directed update [4(| can 
be used both in the grand canonical or canonical modes. 
We describe it in more detail than DQMC since it has 
been developed much more recently. First, the Hamilto- 
nian is written as H = V — T, where V is diagonal in 
the occupation number basis of the direct space, and j 
is the remaining non-diagonal part. The algorithm sam- 
ples an extended partition function represented in the 
interaction picture, 



X t(T)dT 



(2) 



where T T is the time-ordering operator and the time de- 
pendence of any operator A(j) is given by the interaction 
representation A(t) = e rV Ae~ rV . The "Green operator" 
Q is defined by its matrix elements, \ Q \ ^r) = gpq, 
where p and q are the number of creations and destruc- 
tions of particles that transform the (occupation number) 
state ipn into tpL- The matrix g pq is a decreasing func- 
tion of (p + q) [39]. The extended partition function J2]) 
is expanded in T yielding operator strings of the form: 

• ■■t(T L+2 )f(T L + 1 )t(T L )g(T)t(T R )t(T R „ 1 )t(T R „ 2 ) ■ ■ ■ 

(3) 

The Green operator updates the operator string ([3]) by 
propagating in imaginary time, while inserting and re- 
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moving T operators. From its definition, the Green op- 
erator contains terms that do not conserve the number 
of particles. However if the Hamiltonian commutes with 
the operator that measures the total number of particles, 
J\, for example the Hamiltonian ([1]), then only conserva- 
tive terms of Q have non vanishing contributions since the 
trace imposes the same number of particles both at the 
begining and the end of the operator string. As a conse- 
quence the number of particles remains strictly constant 
and the SGF algorithm works in the canonical ensem- 
ble by nature. However, a simple trick can be used in 
order to simulate exactly the grand-canonical ensemble. 
The idea is to add a non conservative part W nc to the 
Hamiltonian, 
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where 7 is an optimization parameter, and allows the 
Green operator to insert at most one H nc operator in the 
string. This allows the number of particles to fluctuate, 
while the addition of the usual term —fjj\ to the Hamil- 
tonian determines the mean number of particles via the 
chemical potential /1. When measuring physical quanti- 
ties, ignoring configurations in which the operator string 
contains a T-L nc operator, corresponds to integrating over 
these configurations the probability of going from a given 
configuration with N particles and no T-Lnc to another one 
with M particles and no "H nc , via intermediate configura- 
tions with 'Hnc. This integrated probability corresponds 
exactly to the probability of going from one configuration 
with N particles to another one with M particles. As a 
result, the configurations of the grand-canonical partition 
function Tr ^-pCM—p 1 *) are generated with the correct 
Boltzmann weight. 

To use this algorithm to simulate fermions in one di- 
mension, we first use the Jordan- Wigncr transformation 
to map the system onto a system of hard core bosons. 
Consequently, this algorithm does not suffer from the 
sign problem but, clearly, it cannot be used in higher di- 
mensions where the Jordan- Wigner transformation fails 
to solve the sign problem. We used this algorithm in the 
canonical ensemble mainly for the confined system where 
it is very convenient to control the populations directly 
rather than tune the chemical potentials. A typical sim- 
ulation for L = 120 and /3 = 32 runs for about 5 days on 
a desktop computer and yields an error of the order of 
0.5% for the peak of the pair momentum distribution.. 

As part of our code verification, the grand canonical 
SGF and DQMC were compared for the same parame- 
ters and found to yield the same results within the error 
bars. Furthermore, we verified that DQMC and grand 
canonical SGF both agreed with the canonical SGF when 
the chemical potentials in the grand canonical cases were 
tuned to give the same populations as the canonical cases, 
see for example Ref.|2g|. The results presented below 
were obtained with these three algorithms; the choice be- 
ing dictated by the specific measurement we were after. 
Because of the equivalence of the three algorithms, we 



will not comment below on which results were obtained 
with which algorithm. 

Using these algorithms, we calculate both the real 
space Green functions of the species, G a , and the pair 
green function, G pa i r , 



G*{1) = Q +la c jtr ) 
G pair (0 = (At +J Aj), 

A 3 = c j'2Cjl, 



(5) 

(6) 
(7) 



where Aj destroys a pair on site j. The Fourier transform 
of G a (l) yields the momentum distributions n a (k) and 
the transform of G pa ; r (7) leads to the pair momentum 
distribution, n pa i r (fc), a central quantity in this work. 
(4) In the non-interacting limit, the Fermi momentum of a 
population is given by kp a = % X T"> wnere L * s the 
length of the system and N a the number of particles. 



III. PHASE DIAGRAM OF THE 
HOMOGENOUS SYSTEM 

We begin with a discussion of the physics in the ab- 
sence of a confining potential, Vt — 0. 

As discussed in the introduction, it is now generally 
agreed that the ground state of a Fermi system with at- 
tractive interactions and imbalanced populations is the 
FFLO state. The mismatch in the Fermi momenta re- 
sults in pair formation with nonzero center-of-mass mo- 
mentum k = ±\kp 1 — kp 2 \. Consequently the pair mo- 
mentum distribution, the Fourier transform of the pair 
Green function Eq.Q, peaks at this value of the mo- 
mentum. This peak at nonzero momentum serves as the 
principal dia gno stic indicating the presence of the FFLO 
state [26|, [28t43l| . 

The situation at finite temperature, which is impor- 
tant experimentally, is less clear. Approximate methods, 
such as mean field, do not always yield the same phase 
diagram. In this section we will map out the phase di- 
agram in the polarization-temperature plane where the 
polarization is defined by, 



P = 



N 1 +N 2 : 



(8) 



where N% (N2) is the majority (minority) population and 
N = Ni +N2 is the total number of particles. To this end, 
we study, at fixed P, the behaviour of the pair momentum 
distribution, rt pa ; r (fc), as a function of the temperature T. 
For very low T, n pa j r (fc) peaks at k ^ and the system is 
in the FFLO state. As T is increased, the peak in rc pa ; r 
gets lower and shifts to k = 0. The temperature at which 
this first happens is the cross-over temperature, T c . Note 
that in this one-dimensional system, transitions at finite 
temperature are not true phase transitions, but rather 
cross-overs. 

This is illustrated in Fig.[T]where we show QMC results 
for n pa i r (fc) as a function of k for several values of the 
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inverse temperature (3. The simulations were done for 
fixed populations, N\ = 13 and N 2 = 7 on a system 
with L = 32 lattice sites and an attractive interaction 
U = — 3.5i. The figure shows clearly that as (3 decreases 
from j3 = 32, the height of the FFLO peak decreases 
and, in fact, shifts to lower k values. The shift to lower k 
values is made more evident by simulating larger systems 
since this gives more k grid points. When the peak at 
nonzero k is equal to n pa i r (0) to within 1%, we consider 
the peak to have shifted to k = and the FFLO state to 
have disappeared. In the Fig. Q] this happens for f3 c « 6. 
Reducing /3 further leads to continued decrease of the 
height of the peak, which remains at k = 0. 

The question then arises as to how the FFLO peak be- 
haves as a function of U at fixed /3, N± and N 2 (and con- 
sequently fixed P). Clearly, for U = there is no pairing 
and no FFLO peak; then, as \U\ is increased, the peak 
at nonzero k forms and its height increases. However, 
we found that as \U\ continues to increase, the FFLO 
peak will reach a maximum height and then start to de- 
crease. We also found that the peak is more sensitive to 
\U\ at higher T. We believe the reason the peak starts 
to go down at high values of \U\ is that with increasing 
attraction, the pairing becomes increasingly localized in 
space and eventually the paired fermions form a very 
tightly bound bosonic molecule and the system resem- 
bles closely a usual Bose-Fermi mixture which does not 
exhibit FFLO peaks. 




FIG. 1: (color online) The effect of temperature on the pair 
momentum distribution. A peak at non-zero momentum is 
a signature of the FFLO state. As j3 decreases, the peak 
disappears at a crossover value f3 c — 1/T C . In this case, j3 c w 
6. The error bars are of the order of the symbol size. 

The peak of n pli i r (k) at fc p0 ak =/= means that the sys- 
tem is, in fact, not homogenous: The spatial pair Green 
function, Eq.©, oscillates as a function of distance with 
a period given by 2tt/ | fc pea k | ■ These oscillations have been 
discussed, for example, in the context of mean field the- 
ory [23|, |4l|. Physically, they indicate that the system 
has regions which are rich in pairs separated by regions 
poor in pairs but rich in the excess fermion species. Such 



oscillations are shown in Fig. [2] for three values of (3. It 
is seen that as the temperature increases and the height 
of the FFLO peak decreases, the oscillations decrease in 
amplitude and eventually disappear as homogeneity is 
restored in the system. 
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FIG. 2: (color online) Pair Green function calculated at low 
temperatures where the system is in the FFLO state (/3 = 32 
and j3 = 16) and at the temperature at which FFLO disap- 
pears (f3 = 6). Note that the oscillations disappear at the 
higher temperature. 

The question then arises as to the nature of the phase 
at T > T c . Two possibilities are: (1) At T > T c the pairs 
are broken and the system is a mixture of two Fermi liq- 
uids or (2) pairs are still present but the system has been 
homogenized by thermal agitation. A first indication is 
given by the energy scales involved. The binding energy 
of the pairs at very low temperature is of the order of 
\U\ and in our system here \U\ = 3.5t. So, to break the 
pairs, an equivalent amount of thermal energy is needed 
which means (3 ss 1/|E/|. For the case discussed in Fig. [5J 
the crossover from FFLO to the uniform phase happens 
at (3 C sa 6 not t/\U\ = 0.286. This means that T c is more 
than an order of magnitude smaller than the temperature 
needed to break the pairs. This, then, favors the conclu- 
sion that when FFLO first disappears, the pairs have not 
yet been broken and the system is in a homogeneous po- 
larized paired phase (PPP). Another piece of evidence is 
provided by studying the average double occupancy of 
the sites given by 

D=(n tl n l2 ) = (AlA t ). (9) 

In the absence of pairing, (nnrii 2 ) = (nn) (rii 2 ) = 
N1N2/L 2 while if pairing is perfect, i.e. if all the mi- 
nority particles are paired, (niin^) = N2/L. We define 
the normalized double occupancy by 

_ D — n\ni , 

V= — , (10) 

n 2 - nin 2 

where n\ — N\/L and n 2 = N 2 /L and we recall that 
A^2 < N\. With this normalization, we have < T> < 1. 
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This quantity is shown in Fig. [3] for three polarizations. 
One can see that the pairing drops significantly at rather 
high temperatures, f3 « 1/U. Thus we conclude that the 
pairs are not broken when the FFLO peak disappears 
but the systems is in a PPP. As one can see, this pairing 
parameter does not saturate for the case we presented. 
One could expect it to reach the maximum value at a very 
strong U limit in the low temperature regime, where both 
thermal and quantum fluctuations are absent. 
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FIG. 3: (color online) Double occupancy versus j3 for three 
polarizations exhibits a very sharp drop for /3 < 1 indicating 
that pairs are being broken near /3 ~ l/|f/|. 

We are now ready to apply the above considerations to 
determine the phase diagram in the (P, T) plane. To this 
end, we keep the total population, N, constant and for 
different values of the polarization, P, determine the tem- 
perature, T c , at which the peak in n pa i r (fc) shifts to k = 0. 
Figure 0] shows the resulting phase diagrams for N = L 
(half filling) and N = L/2 (quarter filling). We mention 
again that the phase boundaries represent crossover be- 
haviour, not phase transitions since this one-dimensional 
quantum system does not have phase transitions at fi- 
nite temperature. The Fermi temperatures shown in the 
figure arc calculated assuming equal populations using 
cf = tkp where ep is the Fermi energy. 

The phase diagrams show clearly that the FFLO phase 
is quite robust, persisting over a wide range of polar- 
izations and to rather high temperatures. For N = L, 
it persists up to T/T F 0.2 and for N = L/2 up 
to T/Tp ~ 0.8. We also see that, in both cases, the 
crossover temperature increases with the polarization up 
to a maximum value after which it decreases again. This 
can be understood physically as follows: When P is 
small, the Fermi "surfaces" of the two populations are 
so close to matching that very little thermal energy is 
needed to get them to match. Thus even at very low 
finite temperature, pairing takes place at zero center-of- 
mass momentum. Therefore, larger polarizations have a 
stabilizing effect on the FFLO phase. 

There arc numerical difficulties with the determination 
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FIG. 4: (color online). Phase diagram in the polarization- 
temperature plane, (a) the system at N = L (half filling), (b) 
the system at N = L/2 (quarter filling). At P=0 system is 
in the BCS state. The regions above the open (red) squares 
and to the left of the vertical (blue) lines are unexplored. 
Up to the open squares the system is in the FFLO phase. 
Phase boundaries represent cross-over behaviour not phase 
transitions. In the FFLO phase, n pa i r (fc) peaks at k 7^ 0, in 
the PPP phase, the peak is at k = 0. 



of the phase diagram at very low and very high polariza- 
tions. At very high polarization, there is a very small 
number of minority particles, making the FFLO signal 
difficult to discern clearly. We were therefore not able 
to examine P > 0.9 in our simulations. For that reason, 
some symbols on the phase boundary at high polariza- 
tion are open, indicating that up to this polarization, the 
system is in the FFLO phase. The solid symbols demark 
the true boundary between FFLO and PPP. 

In addition, at small polarization, very low temper- 
ature is needed to observe the FFLO phase. However, 
there is a practical limit on how low a temperature we 
can simulate since the lower the temperature the longer 
the simulation time needed to obtain precise results. The 
dashed (green) line connecting the origin to the first nu- 
merical points simply schematizes the expected position 
of the boundary. 

A phase diagram in the (P, T) plane was calculated in 
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Ref. [22| using mean field theory (MFT). The general 
shape of the FFLO phase obtained there (Fig. 9 of Rcf. 
[22I ] ) is similar to what we found here. However there are 
very important differences. For example, unlike MFT, 
we have found that there is no direct transition from the 
FFLO phase to the Fermi liquid phase (broken pairs): 
The FFLO phase is destroyed at a temperature which is 
much lower than that required to break the pairs and is 
replaced by the PPP. 

Reference [22[ predicts phase separation between the 
FFLO and PPP at low P and T (Fig. 9 in Ref. 0). 
In order to examine this possibility, we study the density 
histograms in the grand canonical ensemble. The idea 
is as follows: Starting in, say, the PPP, we increase the 
polarization by tuning the chemical potentials, /ii and 
/i2- For each choice of \i\ and /12, we accumulate the 
histograms of the particle populations. If phase separa- 
tion is present, then as the system approaches the phase 
separation region, the density histogram of each species 
should develop a double peak structure. If no such struc- 
ture develops, it means that there is no phase separation 
as the system crosses from the PPP to the FFLO. We first 
verify the correct behaviour of the obtained histogram as 
the size of the system is changed. In Fig. [5] we show the 
histograms for two system sizes at half filling but with 
all other parameters fixed. We see that the histograms 
for the two system sizes agree very well; the main differ- 
ence is that the larger system size (obviously) allows for 
a finer grid of densities which redistributes the values a 
little and exhibits the main peak more clearly. In Fig. [5] 
we show, for fixed inverse temperature (3 = 16, the his- 
tograms for three cases at half filling: In the top panel 
the system is just inside the PPP phase, the middle panel 
the system is at the PPP-FFLO boundary and the bot- 
tom panel the system is just inside the FFLO phase. No 
double peak structure develops, which leads us to con- 
clude that there is no phase separation. This was done 
for several temperatures at low polarization. 

It is useful here to comment on the algorithm choice for 
calculating the histograms. Although the DQMC algo- 
rithm is grand canonical and thus allows for particle num- 
ber fluctuations, it is not useful for calculating the density 
histograms. The reason is that in DQMC one changes the 
realization of the auxiliary Hubbard-Stratonovich field; 
but for each such realization, the fermions have been 
traced over all their possible configurations. On the con- 
trary, in the grand canonical version of the SGF algo- 
rithm, the update is done over the fermion configurations 
themselves. So, the particle number can be measured 
configuration by configuration. 



IV. TRAPPED SYSTEM AT FINITE 
TEMPERATURE 

Continuing earlier work in higher dimension [ToL [Tl| , 
the Rice group 12[ recently reported on experiments in 
one dimensional confined Fermi systems ( 6 Li atoms) with 
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FIG. 5: (color online). Majority (Ni) and minority (N2) den- 
sity histograms at the FFLO-PPP boundary for two system 
sizes. The larger system size offers more grid points and, 
therefore, a finer resolution of the density fluctuations. A sin- 
gle peak is seen for each population indicating the absence of 
phase separation. 
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FIG. 6: (color online). Histograms of local densities for 
L = 30, f3 — 16t and U = — 3.5t. Left peaks correspond 
to the minority species and right ones to the majority. The 
peaks move smoothly as pi and /i2 are tuned to take the sys- 
tem across the boundary between PPP and FFLO. No double 
peak structure is observed, indicating the absence of phase 
separation. 



imbalanced populations. These experiments were done 
in the continuum, i.e. without an optical lattice, and 
focused on the behaviour of the system in three polar- 
ization regimes by measuring the density profiles of the 
fcrmionic species. It was found that the central part of 
the system is always partially polarized whereas the be- 
haviour of the outlying regions depends on the total po- 
larization. For low polarization, P = 0.05, the outlying 
regions were found to be fully paired in the sense that 
the density profiles of the two fermion species matched 
within experimental precision. For medium polarization, 
P = 0.15, the density profiles indicated that the whole 
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system is partially polarized. Finally, for large polariza- 
tion, P = 0.59, the wings were found to be populated 
exclusively by the excess fermion species and thus were 
fully polarized. 

The experiment fl2l ] consisted of a two-dimensional 
array of elongated (one-dimensional) tubes. Along the 
tube, the axial direction, the atoms were confined with 
a trap frequency uj z = 2n x 200Hz; in the central tube, 
the total number of atoms at zero polarization was ap- 
proximately 250 and the temperature was estimated at 
T/Tp rj 0.1. The pair binding energy, e = f? jma\p, 
(where a\p, is the effective one-dimensional scattering 
length), was estimated to be e/e F ~ 5.3 with the Fermi 
energy calculated assuming a balanced system with a to- 
tal of 250 particles. 

In this section, we present QMC results for the 
fermionic Hubbard model, Eq. (JTJ) , in the presence of the 
confining trap. Our goal is to make contact with the 
above mentioned experiment in the continuum; to this 
end we simulate lattice systems that are dilute enough 
so that the fermions in the center of the trap are far 
from forming a flat plateau corresponding to a band in- 
sulator. We introduce the trapping potential in Eq.([T]) 
V T = 0.0007i which corresponds to hw z = The 
total number of particles in our simulations for balanced 
populations is 78, to be compared with 250 in the ex- 
periment. We performed our simulations in the temper- 
ature range 0.008 < T/T F < 0.25 which includes the 
temperature at which the experiments were performed, 
T/Tp = 0.1. In addition, to place our system in the 
same coupling parameter regime as the experiments, we 
present our results for a coupling strength of U = —lOt. 
U is the "pair binding energy" and the value we have cho- 
sen gives \U\/ep = 4.8, close to the experimental value. 

First we look at the system at low temperature when 
we increase the polarization. As in the homogenous case, 
the pair momentum distribution exhibits a maximum at 
fcpeak 7^ (Fig. [7]) and as before fc pea k increases with 
growing polarization. Looking at the density profiles one 
immediately observes that the low and high polarization 
regimes differ significantly. 



A. Low polarization 

In Fig. [5] we show the density profiles at very low tem- 
perature, T/Tp = 0.008, for a system at very small po- 
larization, P = 0.05, corresponding to the open (red) 
squares in Fig. [7] The central region of the system 
is clearly partially polarized: the profiles do not over- 
lap. This partial polarization, i.e. population imbal- 
ance, causes FFLO pairing to take place as is evidenced 
by the pair momentum distribution in Fig. [7l However, 
in a very narrow interval at the edges of the system, the 
density profiles match very closely and the system is fully 
paired [12], 13211 . This fully paired region was studied in the 
continuum [27| using Path Integral Monte Carlo (PIMC) 
simulations. It was shown in Ref. 2jJ that at a po- 
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FIG. 7: (color online). Pair momentum distribution in a 
trapped system for several polarizations. The FFLO peak 
moves to higher momentum values as P increases as in the 
uniform system. The majority population is fixed Ni = 39, 
T/T F = 0.008 in the balanced case Ni = N 2 = 39. 
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FIG. 8: (color online). Density profiles of the two species for 
the low temperature and low polarization case. Very narrow 
fully paired regions are seen at the edges of the cloud. N\ — 39 
and N 2 = 35. 



larization P = 0.04, the fully paired region exists for 
T < 0.025T F and T < 0.035T F for the two strong cou- 
plings studied (the first value corresponds to the smaller 
coupling). The narrowness of this region on the lattice 
was studied at T = in Ref. [H SI]. 

The difference in the density profiles, Fig.O of the two 
species shows regular oscillations indicating that the par- 
tially polarized region is not uniform. Such oscillations 
have been seen before [2(| [28|, [3l| . They correspond to 
the standing wave of length A = 27r/fc pca k, as can be eas- 
ily verified from Figs. [7] and [9l and describe the length 
scale at which the system passes from pair-rich to pair- 
poor regions. This is a striking visual demonstration that 
the FFLO phase is not uniform. 

Next, we examine temperature effects on the system 
at low polarization. The system at /3 — 64, Fig. [51 
is now heated to /3 = 28 (T/T F = 0.016) and /3 = 5 
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FIG. 9: (color online). Density profile differences for different 
polarizations. The oscillations are standing waves showing 
the pair rich (minima) and pair depleted (maxima) regions 
in the confined system. T/Tf — 0.008 for the balanced case 
with N = 78 particles. 
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FIG. 10: (color online). Trapped system at finite temper- 
ature and low polarization, (a) and (b) density profiles at 
T/T F = 0.016 and T/T F = 0.1 respectively, (c) density pro- 
file difference for different temperatures. 



(T/Tp = 0.1) as shown in Fig. [TU] (a) and (b) respec- 
tively. As the temperature increases, the clouds spread 
out and the profiles become more rounded. Figure [101 (c) 
shows what happens to the population differences as the 
temperature rises. As T increases, the population differ- 
ence vanishes in the wings more gradually than at the 
lower temperatures. In addition, the oscillations which 
indicate the presence of FFLO get smoothed out sub- 
stantially at P = 28 and have essentially disappeared 
for /8 = 5. This is confirmed by the behaviour of the 



pair momentum distribution, displayed in Fig. 111! which 
shows the FFLO peak disappearing as T is increased to 
T = 0.0167V. This value is smaller than, but consis- 
tent with, the phase diagram in Ref. [H| (Fig- 1) which 
shows that at P = 0.05 the FFLO phase disappears for 
T > 0.057V. Our value is also consistent with that found 
in Ref. 



27] 



This illustrates the fragility of the fully 
paired paired wings and the FFLO phase at low polar- 
ization. 
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FIG. 11: (color online). Pair momentum distribution of the 
system at P=0.05 (Fig. I10[) and increasing temperature. The 
FFLO peak disappears at T/T F = 0.016. 



B. High polarization 



In section III, we calculated the phase diagram of the 
uniform system and showed that the FFLO phase persists 
to higher temperatures for larger polarization. We now 
demonstrate the same effect in the trapped system. 

Figure [T2l shows the density profiles in the system with 
P = 0.56 at /3 = 64 (T/T F = 0.01). It is clear at this 
large P that the central region is partially polarized and 
the wings are fully polarized, populated only by the ma- 
jority species as was also observed in [28j . The population 
difference in the central region of the system in Fig.[T2"lis 
almost constant but with oscillations, A = 27r/fc pea k, due 
to the presence of FFLO pairing as shown clearly by the 
pair momentum distribution in Fig. [T3] (c). 

As the temperature is increased, this shell structure, 
i.e. partially polarized core and fully polarized wings, 
persists, as can be seen in Fig. [15] (a) and (b). However, 
we observe that the partially polarized core expands in 
size and the fully polarized population in the wings de- 
creases. We also see that the FFLO phase at this large 
polarization is stabilized significantly compared with the 
P = 0.05 case (Fig. [TT]). For P = 0.56, the FFLO peak 
persists, albeit weakly, up to T/Tp = 0.11, vanishing 
completely at T/Tp = 0.25. This result is encouraging 
for the experiments [l2[ which can be done at high polar- 
ization and T/Tp ss 0.1. However, at this temperature, 
the FFLO peak is not very pronounced and might still 
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FIG. 12: (color online). Trapped system at low temperature 
and high polarization. Ni = 39, N 2 = 11, L = 120, U/t = 
— 10, Vt = 0.0007. Tf is calculated for the total balanced 
population of N = 50. Large fully polarized regions are seen 
at the boundaries of the system. 



be difficult to observe experimentally. 

We note that, as in the case of the uniform system, 
the FFLO phase disappears at a temperature which is 
much lower than the contact potential energy. In the 
case above, the FFLO phase disappears at fit sa 3 while 
the contact interactions is U = — lOt. The approximate 
condition to break the pairs is j3U ps 1; we see that /3t = 
3 is not a high enough temperature to break the pairs. 
Therefore, as T is increased, the FFLO phase is replaced 
by the PPP discussed in the previous section and not 
by the normal state composed of the two unpaired spin 
populations. Our QMC result is in disagreement with 
mean field predictions that as T is increased for high 
polarization, the system passes directly from the FFLO 
phase to the normal state [32[ . 

As in the uniform case, the question arises as to 
whether one can stabilize the FFLO phase at higher tem- 
perature simply by increasing the attractive interactcion. 
The answer is the same as in the uniform case: As the 
attractive interaction is increased, the FFLO peak first 
increases in height but then saturates and starts to de- 
crease. In addition, as \U\ is increased, the partially po- 
larized core region shrinks and the polarization in that 
zone increases. This has the effect of shifting the FFLO 
peak to larger values of k. For the fillings we conisdered 
here, the value of the interactions we took, U = — lOt, 
appears to be near optimal. 



V. CONCLUSIONS 

In this paper, we used three QMC algorithms (DQMC, 
canonical SGF and grand canonical SGF) to study the 
behaviour of one-dimensional unbalanced fermion sys- 
tems with attractive interactions governed by the Hub- 
bard Hamiltonian. We explored both the uniform and 
the confined cases. 

In the uniform situation, we mapped the phase dia- 
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FIG. 13: (color online). Same trapped system as in Fig. [12] 
at higher temperature. P — 0.56, L = 120, U/t — —10, 
V t = 0.0007t, Ni = 39, N 2 = 11. (a) and (b) show the 
density profiles and differences for two temperature, (c) Pair 
momentum distribution. 



gram in the polarization-temperature plane at two val- 
ues of the total density (half filling and quarter filling) 
at the same interaction strength U = — 3.5<. Both cases 
show the same general features: (1) The FFLO phase 
is very robust in the ground state and exists for a very 
wide range of polarizations; (2) at small polarization, the 
FFLO phase is very sensitive to temperature increase; 
(3) the FFLO phase persists to higher temperature when 
the polarization is larger and (4) the temperature at 
which the FFLO phase disappears is not high enough 
to break up the pairs leading to a spatially uniform po- 
larized paired phase. 

In the confined case, we performed our simulations 
for system parameters comparable to those in the ex- 
periment of Ref. [HI]. We found the stability features 
of FFLO in the confined phase to be similar to those 
in the uniform case. At low polarization, the FFLO 
phase is destroyed even for a temperature as low as 
T/Tp = 0.016. However, and this is significant for exper- 
imental efforts to detect FFLO in the confined system, we 
found that at large polarization, the FFLO phase persists 
at T/Tp > 0.11 a temperature higher than that of the 
experiment, which has T/Tp rs 0.1. Finally, the temper- 
ature at which the FFLO phase is destroyed is not high 
enough to break the pairs in disagreement with mean 
field results 1321. 
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